Articulation de différents langages (R, JavaScript et Python) pour la géovisualisation avec Quarto

SAGEO, 2023 - Québec, Canada

Nicolas Lambert, Timothée Giraud, Matthieu Viry, Ronan Ysebaert

07 Jun 2023

Python

  • Python : un langage polyvalent, interprété et multi-paradigme

  • De plus en plus utilisé pour la data science

  • Un écosystème robuste pour

Écosystème pour le géospatial

Écosystème pour le géospatial - vecteur

  • Bindings Python de GDAL/OGR
  • Fiona
  • Shapely (bindings Python de GEOS)
  • Pyproj (bindings Python de PROJ)
  • Geopandas - Étend les DataFrames de pandas pour

Écosystème pour le géospatial - GeoPandas

“GeoPandas is an open source project to make working with geospatial data in python easier. GeoPandas extends the datatypes used by pandas to allow spatial operations on geometric types. Geometric operations are performed by shapely. Geopandas further depends on fiona for file access and matplotlib for plotting.”

Format des GeoDataFrame

Écosystème pour le géospatial - raster

  • Rasterio

  • Rasterstats

Écosystème pour le géospatial - Rasterio

Écosystème pour le géospatial - autre

  • Binding Python pour GRASS + Intégration dans les notebooks Jupyter

Écosystème pour le géospatial - Analyse spatiale

Ressources Python Géospatial

  • Yes

  • No

  • Maybe

Interaction R / Python dans Quarto

Python, Quarto et interactivité

  • Utilisation des widgets Jupyter possible (seulement is utilisation de l’engin de rendu jupyter, pas avec knitr - i.e. l’inverse des htmlwidgets en R qui ne fonctionne que si knitr est utilisé)

Python et géospatial dans Quarto

  • Les imports nécessaires :
import numpy as np
import rasterio as rio
import pandas as pd
import geopandas as gpd

from rasterio.warp import calculate_default_transform, reproject, Resampling
from matplotlib import pyplot as plt
from rasterio.plot import show
from rasterstats import zonal_stats
  • On ouvre les jeux de données :
reg_fp = './geom/regio_l.shp'
mun_fp = './geom/munic_s.shp'
lc_fp = '../../Téléchargements/T01_PROVINCE.tif'
lc_categories_fp = '../../Téléchargements/correspondance_raster_CL_COUV.dbf'
dst_epsg_code = 6623

regio_s = gpd.read_file(reg_fp).to_crs(epsg=dst_epsg_code)
munic_s = gpd.read_file(mun_fp).to_crs(epsg=dst_epsg_code)

Python et géospatial dans Quarto

  • Séléction de la zone qu’on souhaite étudier, en utilisant une variable de l’environnement R :
# sel = r.sel
sel = 'Québec'
  • On extrait les entités correspondantes et on les affiche, relativement au reste de la province :
extract = munic_s[munic_s.MUS_NM_MUN == sel]

ax = munic_s.plot(color="lightblue", edgecolor="grey", alpha=0.5)
ax = regio_s.plot(ax=ax, color="orange", alpha=0.8)
ax = extract.plot(ax=ax, color="red")
ax

Python et géospatial dans Quarto

  • Ouverture et reprojection d’un jeu de données raster
dst_crs = f'EPSG:{dst_epsg_code}'

with rio.open(lc_fp) as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
    })
    
    # Read first band
    band = src.read(1)

    # Create the destination array
    img = np.zeros_like(band)

    # Reproject and store the result in the
    # destination array
    reproject(
        band,
        img,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=transform,
        dst_crs=dst_crs,
        resampling=Resampling.nearest,
    )

    # Store the extent of our raster data (we will need it later for plotting)
    left = transform[2]
    top = transform[5]
    right = left + transform[0] * width
    bottom = top + transform[4] * height
    extent = [left, right, bottom, top]
(array([[255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       ...,
       [255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255]], dtype=uint8), Affine(50.00000000000006, 0.0, -830584.8693000015,
       0.0, -50.00000000000006, 1303502.274700009))

Python et géospatial dans Quarto

# Replace '255' values by '0'
img[img == 255] = 0

# Create an empty figure
fig, ax = plt.subplots(figsize=(16, 12))
# Plot the raster using the previously computed extent
show(img, ax=ax, extent=extent, cmap="Set3")
# Plot Québec municipality on top
extract.plot(ax=ax, color='red', edgecolor='red', linewidth=2)

Python et géospatial dans Quarto

# Open the DBF file that contains the correspondences between the codes and the names of land use categories:
categories = pd.DataFrame(gpd.read_file(lc_categories_fp, encoding='utf-8')[['ID', 'CL_COUV', 'Descriptio']])
categories.head(11)
      ID CL_COUV                             Descriptio
0    1.0  010000                 Surfaces artificielles
1    2.0  020000                       Terres agricoles
2    3.0  060000             Milieux humides forestiers
3    4.0  070000  Milieux humides herbacés ou arbustifs
4    5.0  100000        Plans et cours d’eau intérieure
5    6.0  110101    Forêts de conifères à couvert fermé
6    7.0  110102     Forêts de feuillus à couvert fermé
7    8.0  110103          Forêts mixtes à couvert fermé
8    9.0  110200                Forêts à couvert ouvert
9   10.0  000000                         Pas de données
10  11.0  000001               En attente de traitement

Python et géospatial dans Quarto

  • Compute the zonal statistics (how many cell of each type in the selection):
stats = zonal_stats(extract, img, affine=transform, categorical=True, nodata=0)
stats
[{1: 85316, 2: 12475, 3: 2737, 4: 756, 5: 3494, 6: 3556, 7: 38023, 8: 38147, 9: 69, 10: 43}, {8: 4}]
  • Lets group the values into a single dict:
# We convert the dict objects into Counter objects, which allows us to add them up by simply using the "+" operator.
from collections import Counter
stats = Counter(stats[0]) + Counter(stats[1])
stats = dict(stats) # Convert back to dict
stats
{1: 85316, 2: 12475, 3: 2737, 4: 756, 5: 3494, 6: 3556, 7: 38023, 8: 38151, 9: 69, 10: 43}

Python

  • Lets plot an histogram of this result, using the appropriate names of land use categories:
# Use the category names from the `categories` DataFrame
x = [
  categories[categories.ID == int(v)]['Descriptio'].iloc[0]
  for v in stats.keys()
]
# and the values we just computed with the `zonal_stats` function
height = stats.values()

fig, ax = plt.subplots(figsize=(8, 5))
bar_container = ax.bar(x, height)
ax.set_xticklabels(
    ax.get_xticklabels(),
    rotation=45,
    horizontalalignment='right'
)
ax

Python